function matrix_square_root()
% MATRIX_SQUARE_ROOT 矩阵平方根计算
% 
% 展示计算矩阵平方根的各种方法和应用

fprintf('=== 矩阵平方根计算 ===\n\n');

%% 1. 正定矩阵的平方根
fprintf('1. 正定矩阵的平方根\n');
fprintf('目标：计算正定矩阵的主平方根\n\n');

% 正定矩阵
A1 = [4, 2; 2, 3];
fprintf('正定矩阵 A1 =\n');
disp(A1);

% 验证正定性
eigenvals1 = eig(A1);
fprintf('特征值: [%.4f, %.4f]\n', eigenvals1);
fprintf('正定性: %s\n', mat2str(all(eigenvals1 > 0)));

% 方法1：特征值分解
[V1, D1] = eig(A1);
sqrt_D1 = diag(sqrt(diag(D1)));
A1_sqrt_eigen = V1 * sqrt_D1 * V1';

fprintf('方法1：特征值分解\n');
fprintf('A1^(1/2) =\n');
disp(A1_sqrt_eigen);

% 验证：A1_sqrt * A1_sqrt = A1
verification1 = A1_sqrt_eigen * A1_sqrt_eigen;
fprintf('验证 (A^(1/2))^2 = A，误差: %.2e\n', norm(verification1 - A1, 'fro'));

% 方法2：Cholesky分解
L1 = chol(A1, 'lower');
A1_sqrt_chol = L1;

fprintf('方法2：Cholesky分解\n');
fprintf('A1^(1/2) (Cholesky) =\n');
disp(A1_sqrt_chol);

% 验证
verification1_chol = A1_sqrt_chol * A1_sqrt_chol';
fprintf('验证 L*L'' = A，误差: %.2e\n', norm(verification1_chol - A1, 'fro'));

% 比较两种方法
fprintf('两种方法的差异: %.2e\n', norm(A1_sqrt_eigen - A1_sqrt_chol, 'fro'));
fprintf('注意：Cholesky给出下三角平方根，特征值分解给出对称平方根\n');

% 方法3：MATLAB内置函数
A1_sqrt_matlab = sqrtm(A1);
fprintf('方法3：MATLAB sqrtm函数\n');
fprintf('A1^(1/2) (sqrtm) =\n');
disp(A1_sqrt_matlab);

verification1_matlab = A1_sqrt_matlab * A1_sqrt_matlab;
fprintf('验证误差: %.2e\n', norm(verification1_matlab - A1, 'fro'));

%% 2. 半正定矩阵的平方根
fprintf('\n2. 半正定矩阵的平方根\n');
fprintf('目标：处理奇异的半正定矩阵\n\n');

% 半正定矩阵（奇异）
A2 = [1, 1, 1; 1, 1, 1; 1, 1, 1];
fprintf('半正定矩阵 A2 =\n');
disp(A2);

eigenvals2 = eig(A2);
fprintf('特征值: [%.4f, %.4f, %.4f]\n', eigenvals2);
fprintf('秩: %d\n', rank(A2));

% 特征值分解方法
[V2, D2] = eig(A2);
eigenvals2_sorted = sort(diag(D2), 'descend');
[~, sort_idx] = sort(diag(D2), 'descend');
V2_sorted = V2(:, sort_idx);

% 处理零特征值
tol = 1e-12;
positive_idx = eigenvals2_sorted > tol;
n_positive = sum(positive_idx);

fprintf('正特征值个数: %d\n', n_positive);

if n_positive > 0
    % 构造平方根
    sqrt_eigenvals = zeros(size(eigenvals2_sorted));
    sqrt_eigenvals(positive_idx) = sqrt(eigenvals2_sorted(positive_idx));
    
    A2_sqrt = V2_sorted * diag(sqrt_eigenvals) * V2_sorted';
    
    fprintf('半正定矩阵的平方根 A2^(1/2) =\n');
    disp(A2_sqrt);
    
    % 验证
    verification2 = A2_sqrt * A2_sqrt;
    fprintf('验证误差: %.2e\n', norm(verification2 - A2, 'fro'));
    
    % 检查平方根的性质
    fprintf('平方根的秩: %d\n', rank(A2_sqrt));
    fprintf('平方根的特征值: [%.4f, %.4f, %.4f]\n', sort(eig(A2_sqrt), 'descend'));
end

%% 3. 矩阵平方根的唯一性
fprintf('\n3. 矩阵平方根的唯一性\n');
fprintf('目标：展示主平方根的唯一性\n\n');

A3 = [5, 2; 2, 2];
fprintf('矩阵 A3 =\n');
disp(A3);

% 主平方根（特征值分解）
[V3, D3] = eig(A3);
A3_sqrt_principal = V3 * diag(sqrt(diag(D3))) * V3';

fprintf('主平方根 =\n');
disp(A3_sqrt_principal);

% 另一个平方根（改变符号）
A3_sqrt_other = V3 * diag([-sqrt(diag(D3, 1)); sqrt(diag(D3, 2))]) * V3';

fprintf('另一个平方根 =\n');
disp(A3_sqrt_other);

% 验证两个都是平方根
verification3_1 = A3_sqrt_principal * A3_sqrt_principal;
verification3_2 = A3_sqrt_other * A3_sqrt_other;

fprintf('主平方根验证误差: %.2e\n', norm(verification3_1 - A3, 'fro'));
fprintf('另一个平方根验证误差: %.2e\n', norm(verification3_2 - A3, 'fro'));

% 主平方根的特征：所有特征值为正
eigenvals3_principal = eig(A3_sqrt_principal);
eigenvals3_other = eig(A3_sqrt_other);

fprintf('主平方根特征值: [%.4f, %.4f]\n', eigenvals3_principal);
fprintf('另一个平方根特征值: [%.4f, %.4f]\n', eigenvals3_other);

%% 4. 矩阵平方根的计算方法比较
fprintf('\n4. 矩阵平方根的计算方法比较\n');
fprintf('目标：比较不同计算方法的性能和精度\n\n');

% 测试矩阵
sizes = [10, 50, 100];
fprintf('矩阵规模    特征值分解    Schur分解    sqrtm函数\n');
fprintf('--------    ----------    ---------    ---------\n');

for n = sizes
    % 生成随机正定矩阵
    A_rand = randn(n, n);
    A_test = A_rand' * A_rand + eye(n);  % 确保正定
    
    % 方法1：特征值分解
    tic;
    [V_test, D_test] = eig(A_test);
    sqrt_test_eigen = V_test * diag(sqrt(diag(D_test))) * V_test';
    time_eigen = toc;
    
    % 方法2：Schur分解
    tic;
    [U_test, T_test] = schur(A_test);
    sqrt_T = sqrtm(T_test);  % 上三角矩阵的平方根
    sqrt_test_schur = U_test * sqrt_T * U_test';
    time_schur = toc;
    
    % 方法3：MATLAB内置函数
    tic;
    sqrt_test_matlab = sqrtm(A_test);
    time_matlab = toc;
    
    fprintf('%6d      %10.6f    %9.6f    %9.6f\n', ...
            n, time_eigen, time_schur, time_matlab);
    
    % 精度比较（仅对小矩阵）
    if n <= 50
        error_eigen = norm(sqrt_test_eigen * sqrt_test_eigen - A_test, 'fro');
        error_schur = norm(sqrt_test_schur * sqrt_test_schur - A_test, 'fro');
        error_matlab = norm(sqrt_test_matlab * sqrt_test_matlab - A_test, 'fro');
        
        fprintf('        精度: %.2e        %.2e        %.2e\n', ...
                error_eigen, error_schur, error_matlab);
    end
end

%% 5. 矩阵平方根的应用：协方差矩阵的白化
fprintf('\n5. 矩阵平方根的应用：协方差矩阵的白化\n');
fprintf('目标：使用矩阵平方根进行数据白化\n\n');

% 生成相关数据
n_samples = 200;
n_vars = 3;
rng(123);

% 原始数据（有相关性）
Sigma_true = [2, 1, 0.5; 1, 1.5, 0.8; 0.5, 0.8, 1];
L_true = chol(Sigma_true, 'lower');
X_original = randn(n_samples, n_vars) * L_true';

fprintf('原始数据维度: %d×%d\n', size(X_original));

% 计算样本协方差矩阵
X_centered = X_original - mean(X_original);
Sigma_sample = (X_centered' * X_centered) / (n_samples - 1);

fprintf('样本协方差矩阵 =\n');
disp(Sigma_sample);

% 白化变换：使用协方差矩阵的逆平方根
[V_sigma, D_sigma] = eig(Sigma_sample);
Sigma_inv_sqrt = V_sigma * diag(1./sqrt(diag(D_sigma))) * V_sigma';

fprintf('协方差矩阵的逆平方根 =\n');
disp(Sigma_inv_sqrt);

% 应用白化变换
X_whitened = X_centered * Sigma_inv_sqrt';

% 验证白化效果
Sigma_whitened = (X_whitened' * X_whitened) / (n_samples - 1);
fprintf('白化后的协方差矩阵 =\n');
disp(Sigma_whitened);

fprintf('白化效果（应该接近单位矩阵）:\n');
fprintf('与单位矩阵的差异: %.2e\n', norm(Sigma_whitened - eye(n_vars), 'fro'));

% 验证逆变换
X_reconstructed = X_whitened * (Sigma_inv_sqrt')^(-1);
reconstruction_error = norm(X_centered - X_reconstructed, 'fro');
fprintf('重构误差: %.2e\n', reconstruction_error);

%% 6. 矩阵平方根的数值稳定性
fprintf('\n6. 矩阵平方根的数值稳定性\n');
fprintf('目标：分析不同条件数下的数值稳定性\n\n');

condition_numbers = [1e2, 1e6, 1e10, 1e12];
fprintf('条件数        计算成功    验证误差      条件数变化\n');
fprintf('--------      --------    --------      ----------\n');

for cond_target = condition_numbers
    % 构造指定条件数的正定矩阵
    n = 4;
    eigenvals_target = [1, 1/sqrt(cond_target), 1/cond_target^(1/3), 1/cond_target];
    V_rand = orth(randn(n, n));
    A_cond = V_rand * diag(eigenvals_target) * V_rand';
    
    try
        % 计算平方根
        A_sqrt_cond = sqrtm(A_cond);
        success = true;
        
        % 验证误差
        verification_error = norm(A_sqrt_cond * A_sqrt_cond - A_cond, 'fro') / norm(A_cond, 'fro');
        
        % 条件数变化
        cond_original = cond(A_cond);
        cond_sqrt = cond(A_sqrt_cond);
        
        fprintf('%.0e      %8s    %.2e      %.2e\n', ...
                cond_target, 'Yes', verification_error, cond_sqrt/sqrt(cond_original));
        
    catch
        fprintf('%.0e      %8s    %8s      %8s\n', ...
                cond_target, 'No', 'N/A', 'N/A');
    end
end

fprintf('\n理论：如果A的条件数为κ，则A^(1/2)的条件数约为√κ\n');

%% 7. 复矩阵的平方根
fprintf('\n7. 复矩阵的平方根\n');
fprintf('目标：处理复数矩阵的平方根\n\n');

% 复Hermitian正定矩阵
A7 = [2, 1+1i; 1-1i, 3];
fprintf('复Hermitian矩阵 A7 =\n');
disp(A7);

% 检查Hermitian性
hermitian_error = norm(A7 - A7', 'fro');
fprintf('Hermitian性检查: %.2e\n', hermitian_error);

% 检查正定性
eigenvals7 = eig(A7);
fprintf('特征值: [%.4f, %.4f]\n', real(eigenvals7));
fprintf('正定性: %s\n', mat2str(all(real(eigenvals7) > 0 & abs(imag(eigenvals7)) < 1e-12)));

if all(real(eigenvals7) > 0)
    % 计算复矩阵的平方根
    A7_sqrt = sqrtm(A7);
    fprintf('复矩阵的平方根 =\n');
    disp(A7_sqrt);
    
    % 验证
    verification7 = A7_sqrt * A7_sqrt;
    fprintf('验证误差: %.2e\n', norm(verification7 - A7, 'fro'));
    
    % 检查平方根的Hermitian性
    sqrt_hermitian_error = norm(A7_sqrt - A7_sqrt', 'fro');
    fprintf('平方根的Hermitian性: %.2e\n', sqrt_hermitian_error);
end

%% 总结
fprintf('\n=== 矩阵平方根总结 ===\n');
fprintf('1. 正定矩阵有唯一的主平方根\n');
fprintf('2. 半正定矩阵的平方根存在但可能不唯一\n');
fprintf('3. 特征值分解是最直接的计算方法\n');
fprintf('4. Cholesky分解给出下三角平方根\n');
fprintf('5. 数值稳定性与矩阵条件数相关\n');
fprintf('6. 在数据白化中有重要应用\n');
fprintf('7. 可以扩展到复Hermitian矩阵\n\n');

fprintf('计算方法选择：\n');
fprintf('• 对称正定矩阵：特征值分解或Cholesky\n');
fprintf('• 一般矩阵：Schur分解\n');
fprintf('• 高精度要求：MATLAB sqrtm函数\n');
fprintf('• 大规模矩阵：迭代方法\n\n');

fprintf('应用领域：\n');
fprintf('• 数据预处理和白化\n');
fprintf('• 协方差矩阵分析\n');
fprintf('• 优化算法\n');
fprintf('• 控制理论\n');
fprintf('• 统计学和机器学习\n\n');

end